#ifndef _RJBUTILITIES_H
#define _RJBUTILITIES_H

#include "itkLabelStatisticsImageFilter.h"
#include <itkBinaryThresholdImageFilter.h>

#if 0
// this uses a modified version of label statistics image filter
/////////////////////////////////////////////////////////
template <class RawIm, class MaskIm>
void
computeMaskQuantiles(typename RawIm::Pointer raw, typename MaskIm::Pointer mask, 
		     float &median,
		     float &q25,
		     float &q75)
{
  // robust rescaling based on median and IQR
  typedef typename itk::LabelStatisticsImageFilter<RawIm, MaskIm> StatsType;
  
  typename StatsType::Pointer statsfilt = StatsType::New();
  statsfilt->SetInput(raw);
  statsfilt->SetLabelInput(mask);
  statsfilt->Update();

  typename StatsType::RealType Mn = statsfilt->GetMinimum(1);
  typename StatsType::RealType Mx = statsfilt->GetMaximum(1);

  statsfilt->SetHistogramParameters(100, Mn, Mx);
  statsfilt->SetUseHistograms(true);
  statsfilt->Modified();
  statsfilt->Update();
  //typename StatsType::RealType median = statsfilt->GetQuantile(1, 0.5);
  //typename StatsType::RealType q25 = statsfilt->GetQuantile(1, 0.25);
  //typename StatsType::RealType q75 = statsfilt->GetQuantile(1,
  //0.75);
  median = (float)statsfilt->GetQuantile(1, 0.5);
  q25 = (float)statsfilt->GetQuantile(1, 0.25);
  q75 = (float)statsfilt->GetQuantile(1,0.75);

}
#endif
template <class RawIm, class MaskIm>
void
computeMaskMeanVar(typename RawIm::Pointer raw, typename MaskIm::Pointer mask,
		   typename MaskIm::PixelType Lab,
		   double &Mean,
		   double &Var)
{
  // robust rescaling based on median and IQR
  typedef typename itk::LabelStatisticsImageFilter<RawIm, MaskIm> StatsType;
  
  typename StatsType::Pointer statsfilt = StatsType::New();
  statsfilt->SetInput(raw);
  statsfilt->SetLabelInput(mask);
  statsfilt->Update();

  Mean = (double)statsfilt->GetMean(Lab);
  Var = (double)statsfilt->GetVariance(Lab);

}

/////////////////////////////////////////////////////////
template <class RawIm, class MaskIm>
typename MaskIm::Pointer doThresh(typename RawIm::Pointer raw, float threshVal, float scale = 1.0)
{
  // a convenience function - sacrifices streaming
  typedef typename itk::BinaryThresholdImageFilter<RawIm, MaskIm> ThreshType;
  typename ThreshType::Pointer wthresh = ThreshType::New();
  wthresh->SetInput(raw);
  // take into account spm's scaling
  wthresh->SetUpperThreshold((typename RawIm::PixelType)(threshVal * scale));
  wthresh->SetLowerThreshold(0);
  wthresh->SetInsideValue(0);
  wthresh->SetOutsideValue(1);
  
  typename MaskIm::Pointer result = wthresh->GetOutput();
  result->Update();
  result->DisconnectPipeline();
  return(result);
}

/////////////////////////////////////////////////////////
template <class RawIm, class MaskIm>
typename MaskIm::Pointer doThresh2(typename RawIm::Pointer raw, float threshVal, float scale = 1.0)
{
  // a convenience function - sacrifices streaming
  typedef typename itk::BinaryThresholdImageFilter<RawIm, MaskIm> ThreshType;
  typename ThreshType::Pointer wthresh = ThreshType::New();
  wthresh->SetInput(raw);
  // take into account spm's scaling
  wthresh->SetUpperThreshold((typename RawIm::PixelType)(threshVal * scale));
  wthresh->SetLowerThreshold(0);
  wthresh->SetInsideValue(1);
  wthresh->SetOutsideValue(0);
  
  typename MaskIm::Pointer result = wthresh->GetOutput();
  result->Update();
  result->DisconnectPipeline();
  return(result);
}

#endif
